Marginalization using the metric of likelihood

Although the likelihood function is normalizeable with respect to the data there is no guarantee that the same holds with respect to the model parameters. This may lead to singularities in the expectation value integral of these parameters, especially if the prior information is not sufficient to take care of finite integral values. However, the problem may be solved by obeying the correct Riemannian metric imposed by the likelihood. This will be demonstrated for the example of the electron temperature evaluation in hydrogen plasmas. INTRODUCTION Given data ~ d, a linear parameter and some function ~ f meant to explain the data, we have ~ d= ~ f(T )+~" : (1) The vectors shall have dimension N according to the number of quantities measured. Due to the measurement process the data is corrupted by noise, where h"i = 0 and h"2i= 2. Then by the principle of Maximum Entropy the likelihood function reads p(Dj ; ; ~ f ;I)/ exp( 1 2 2Xi [di fi℄2) ; (2) which is clearly normalizeable for the data ~ d and bound for every parameter showing up as a functional dependency in f . The situation may change when we are looking for the expectation value of some parameter of f , let say f = f(T ). Then we need to evaluate the posterior of T with hT i / Z T p(dT jD;I) : (3) In order to connect the unknown posterior to the known likelihood we marginalize over all the parameters which enter the problem, that is in our problem and : p(dT jD;I) = Z Z p(dT;d ;d jD;I) ; (4) and make use of Bayes theorem: p(T; ; jD;I)/ p(DjT; ; ;I) p(T; ; jI) : (5) Commonly, the infinitesimal elements in equation (4) are identified with p(dT;d ;d jD;I) = p(T; ; jD;I) dT d d : (6) In mathematical terms this would mean that the probability functions live in euclidean space. They do not. RIEMANNIAN METRIC Parameterizations correspond to choices of coordinate systems. The problem to be solved has to be invariant against reparametrizations [1], i.e. in the space of the probability functions one has to get the same answer no matter what parameters were chosen to describe a model. Therefore one is in need of a length measure which takes care of defining a distance between different elements of this probability function space. This task is done by applying differential geometry to statistical models, an approach which was baptized ’information geometry’ by S. Amari [2]. Eq. (6) then reads correctly p(dT;d ;d jD;I) = p(T; ; jD;I) (dT;d ;d ) : (7) (d~ ) = (~ )d~ is the natural Riemannian metric on a regular model (in our case the model is parameterized by ~ = (T; ; )). It results from second variations of the entropy [2, 3] and is given by (d~ ) =qdetg(~ )d~ (8) where g is the Fisher information matrix: gij = * 2 logp(Dj~ ;I) i j + : (9) For the above likelihood the metric reads explicitly ( ; ;T )/ 3vuuut"Xi f 2 i #24Xi fi T !235 "Xi fi fi T #2 : (10) Notice that this approach is based on the assumption that the hypothesis space of the likelihood defines the metric to be calculated in. This may not be the case if some prior information was already used during data acquisition, e.g. the experimentalist uses his expert knowledge in separating ’correct’ data from the rest. The latter instantly rules out certain parts of all possible realizations of the likelihood function and results in a different hypothesis space. SIMPLE EXAMPLE First we want to demonstrate the relevance of using the correct metric with a simple example which already has all the features of the real world problem further down. fi(T ) = T (T +xi) 1xi ; (11) where the notation in i corresponds to the data points di. For simplification let us assume that the variance 2 is known and we only have to marginalize over in order to get the posterior. What happens if we do not use the Riemannian metric? Then the marginalization integral over reads p(T jD;I)/ Z d p(DjT; ;I) p( jI) : (12) In order to facilitate analytic calculation the exponent of the likelihood is written in a quadratic form over Xi [di fi℄2 = (~ fT ~ f)[ 0℄2+24~ dT ~ d (~ dT ~ f)2 ~ fT ~ f 35 ; (13) where 0 = ~ dT ~ f=~ fT ~ f . For the prior p( jI) the only thing we know is that will be something in between an upper and a lower limit, where it is reasonable to assume that the upper (lower) bound is given by an unknown factor n (1/n) of the value 0 where the maximum of the likelihood occurs. The principle of maximum entropy gives a flat prior with p( jI) = ( 1 n 0 8 0 n 0 0 else : (14) The integral over the c-dependent parts then reads 1 n 0 Z n 0 0 d exp 1 2 2 (~ fT ~ f)[ 0℄2 : (15) One may check that for ~ fT ~ f 2 it is allowed to shift the integral boundaries to += infinity with affecting the value of the integrand up to a small error only. As a matter of fact for the chosen model parameters of N=3, xi=i, T=1, c=1 and =0.1 the error is in the order of 10 7 of the correct integral. Notice that this is almost the same for every T in between 0 and infinity. We finally get p(T jD;I)/ q~ fT ~ f ~ dT ~ f exp8<: 1 2 2 24~ dT ~ d (~ dT ~ f)2 ~ fT ~ f 359=; : (16) A look at the behavior for large and small T gives lim T!0p(T jD;I) / pN Pidi exp( 1 2 2 "~ dT ~ d (Pi di)2 N #) / onst ; (17) lim T!1p(T jD;I) / p~xT~x ~ dT~x exp8<: 1 2 2 24~ dT ~ d (~ dT~x)2 ~xT~x 359=; / onst : (18) 10 −3 10 −2 10 −1 10 0 10 1 10 2 10 3 T 10 −4 10 −3 10 −2 10 −1 10 0